Differential Growth Rates of Northeast US Fishes Following Temperature Regime Shift

Author
Affiliation

Gulf of Maine Research Institute

Published

August 11, 2022

Growth Regimes in the NW Atlantic Following Temperature Regime Shift

Introduction

Literature on TSR and Growth Expectations

Introduce Study Area and Existing Research

Recent changes in Gulf stream positioning have altered the relative influences of the Gulf stream and Labrador current on temperature and salinity regimes for the region. Understanding that this new temperature and salinity regime should alter both the energetics of the region, as well as food-web dynamics. We seek to identify whether these changes have had a measurable impact on the growth of several groundfish species.

Despite many of these species living at or near the bottom, research suggests that the forces driving these thermal regimes should have far-reaching impacts biological impacts. These impacts may directly change the near-bottom environment through mixing patterns, and/or indirectly through food-web impacts.

Through the lens of an environmental regime change we will assess changes to body size, size-at-age, and growth characteristics derived using the von-bertalanffy growth equation.

Methods

Sea Surface Temperature

Global Sea surface temperature data was obtained via NOAA’s optimally interpolated SST analysis, providing daily temperature values at a 0.25° latitude x 0.25° longitude resolution (Reynolds et al. 2007). A daily climatology for every 0.25° pixel in the global data set was created using average daily temperatures spanning the period of 1982-2011. Daily anomalies were then computed as the difference between observed temperatures and the daily climatological average. OISSTv2 data used in these analyses were provided by the NOAA PSL, Boulder, Colorado, USA from their website at https://psl.noaa.gov.

Code
# Load the regional temperatures from {targets}
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(regional_oisst))


# Get regional averages
temp_regimes <- regional_oisst %>% 
  filter(yr > 1981) %>% 
  mutate(
    regime = ifelse(yr < 2010, "1982-2009 Average", "2010-2021 Average"),
    survey_area = ifelse(survey_area == "all", "Full Region", survey_area),
    survey_area = factor(survey_area, levels = c("Full Region", "GoM", "GB", "SNE", "MAB"))) %>% 
  group_by(survey_area, regime) %>% 
  mutate(
    avg_temp_c = round(mean(sst), 2),
    temp_label = str_c(avg_temp_c, "\u00b0C"),
    avg_anom_c = mean(sst_anom),
    #avg_anom_c = ifelse(regime == "1982-2009 Regime", 0, avg_anom_c),
    x = ifelse(regime == "1982-2009 Average", 1981.5, 2009.5),
    xend = ifelse(regime == "1982-2009 Average", 2009.5, 2021.5),
         ) %>% 
  ungroup()

Temperature Study Region

Temperature data was regionally averaged to match the sampling area for the fisheries independent survey program providing the age-at-length data. SST anomalies were averaged over the entire sampling region, consisting of continental shelf habitats from Cape Hatteras to Nova Scotia, to produce a daily time series. This time series was then processed into an annual timeseries of surface temperatures and anomalies. All area-averaging was done with area-weighting to account for differences in the relative areas of the latitude/longitude grid of the OISSTv2 data.

Temperature Regime Shift

Text about identifying the regime shift

Distribution Shift Analyses

Fishery Independent data on groundfish size-at-age was collected as part of the NEFSC’s northeast trawl survey. This survey is conducted each year in the Spring and in the fall, with sample locations determined following a stratified-random survey design with effort allocated in proportion to stratum area. Add an Explanation of the Strata we Exclude. Correction factors were applied for changes in vessels, gear, and doors when appropriate. Prior to sampling, a CTD cast is performed to collect environmental conditions for the water column at/near the start of the trawl sample providing information on bottom habitat conditions like bottom temperatures. Distribution shift analyses were limited to data from the Spring survey season. Previous work has shown no significant changes in timing of sampling, or the mean annual latitude and longitude of the Spring survey sampling (Nye et al. 2009).

What do we calculate and for how many species?

Our analyses followed the continued movements of 30 species consistent with Nye et al. 2009. These species were originally selected for being consistently sampled across all years, and as representatives of a wide range of taxonomic groups.

Code
# Access the tidy survdat Data with {targets}
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(survdat_clean))


# Filter to Just spring?
dist_dat <- survdat_clean %>% 
  filter(season == "Spring")

# Just the Species, or as Stocks?
nye_species <- c(
  "spiny dogfish",
  "little skate",
  "thorny skate",
  "winter skate",
  "alewife",
  "american shad",
  "atlantic herring",
  "atlantic cod",
  "haddock",
  "pollock",
  "silver hake",
  "red hake",
  "spotted hake",
  "white hake",
  "cusk",
  "goosefish",
  "american plaice",
  "atlantic halibut",
  "yellowtail flounder",
  "winter flounder",
  "fourspot flounder",
  "summer flounder",
  "windowpane",
  "atlantic mackerel",
  "ocean pout",
  "atlantic wolffish",
  "blackbelly rosefish",
  "longhorn sculpin",
  "sea raven",
  "acadian redfish"
)


# Filter to those species
dist_dat <- dist_dat %>% filter(comname %in% nye_species)

# Get the total biomass at each station
# Totals up the abundance and biomass across sex*
station_totals <- dist_dat %>% 
  group_by(est_year, survey_area, stratum, tow, est_towdate, avgdepth, bottemp, decdeg_beglat, decdeg_beglon, comname) %>% 
  summarise(
    biomass_kg = sum(biomass_kg, na.rm = T),
    .groups = "drop")

Calculating Spatial Metrics

Movements through time were characterized for each species using the following metrics: center of biomass, area-occupied, mean depth of occurrence, and the mean temperature of occurrence. Each metric was weighted by the biomass of a species sampled at each station such that.

\[X_j = \frac{ \sum_{}^{} w_iX_{ij} }{ \sum~w_i }\]

Where \(X\) is the spatial metric, \(j\) is the year, & \(w\) is the biomass (in kg) for each station \(i\).

Code
# Function to run the group weighted means on different groups using ...
# Does all variables: depth, temp, lat, lon & their variance (sd)
group_spat_metrics <- function(station_totals, ...){
  station_totals %>% 
    group_by(comname, ...) %>% 
    summarise(
      total_biomass   = sum(biomass_kg),
      avg_biomass     = mean(biomass_kg),
      biomass_sd      = sd(biomass_kg),
      avg_depth       = weightedMean(avgdepth, w = biomass_kg, na.rm = T),
      avg_temp        = weightedMean(bottemp, w = biomass_kg, na.rm = T),
      avg_lat         = weightedMean(decdeg_beglat, w = biomass_kg, na.rm = T),
      avg_lon         = weightedMean(decdeg_beglon, w = biomass_kg, na.rm = T),
      depth_sd        = weightedSd(avgdepth, w = biomass_kg, na.rm = T),
      temp_sd         = weightedSd(bottemp, w = biomass_kg, na.rm = T),
      lat_sd          = weightedSd(decdeg_beglat, w = biomass_kg, na.rm = T),
      lon_sd          = weightedSd(decdeg_beglon, w = biomass_kg, na.rm = T),
      .groups = "drop")
}


# Function to calculate weighted mean & median across biomass variables
# returns a common format that can join easily to global benchmarks
biomass_weighted_var <- function(station_totals, var_col, ...){
  var_lab <- deparse(substitute(var_col))
  mean_df <- station_totals %>% 
  group_by(comname, ...) %>% 
  summarise(
    biom_median  = weightedMedian({{var_col}}, biomass_kg, na.rm = T),
    biom_mu      = weightedMean({{var_col}},   biomass_kg, na.rm = T),
    biom_sd      = weightedSd({{var_col}}, biomass_kg, na.rm = T), 
    .groups = "drop") %>% 
    mutate(var = var_lab)
  
  return(mean_df)
  
}



# 1. Yearly Spatial Metrics
year_avgs <- group_spat_metrics(station_totals, est_year) 


# 2. Average Metrics Across Temperature Regimes
regime_avgs <- station_totals %>% 
  filter(est_year >= 2000) %>% 
  mutate(temp_regime = ifelse(est_year < 2010, "Early Regime", "Warm Regime")) %>% 
  group_spat_metrics(., temp_regime) 


# Variance across all years within each species
global_benchmarks <- group_spat_metrics(station_totals) 

Analysis of Spatial Metrics

Size at Age Analyses

The available age data is a subset of the overall catch data, which contains additional information on individual fishes such as otolith-derived aging that require additional workup to ascertain. Not all fish species have age available for the full time period, with special attention being given to aid in the management of commercially valuable species.

Code
# Access Biological Data with {targets}
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(survdat_biological))

# 1. Prepare Bio Data for Age analysis for regimes:

# Drop NA age data, restrict to summer and fall
nmfs_bio <- survdat_biological %>%
  filter(is.na(age) == FALSE,
         season %in% c("Spring", "Fall"))%>% 
      as.data.frame()
  

#Add regime and decade labels
nmfs_bio <- nmfs_bio %>% 
  mutate(
    yearclass = est_year - (age-1),
    regime = ifelse(est_year < 2010, "Early Regime", "Warm Regime"),
    decade = floor_decade(est_year)) 

rm(survdat_biological)



####  Growth Data Prep  #### 

# Filter to just the data for the two regimes
regime_dat <- filter(nmfs_bio, est_year >= 2000, est_year <= 2019)


# Rank species by how many measurements there are
species_abunds <- regime_dat %>% 
  count(comname) %>% 
  arrange(desc(n)) # ordered by number measured

# # Reorder alphabetically And take
# vonbert_species <- sort(species_abunds$comname[1:17])


# Drop the ones that don't work
vonbert_species <- species_abunds %>% 
  filter(comname %not in% c("goosefish", "fourspot flounder", "cusk", "ocean pout", "bluefish")) %>% 
  arrange(comname) %>% 
  pull(comname)

Size at age analyses were performed on 18 species. Species were omitted from analyses if age data was not available across both temperature regimes (bluefish & ocean pout) or if their physiology prevented accurate aging (elasmobranch species). The following 18 species had sufficient age data across both temperature regimes to be included in the analysis: acadian redfish, american plaice, atlantic cod, atlantic herring, atlantic mackerel, black sea bass, butterfish, haddock, pollock, red hake, scup, silver hake, summer flounder, white hake, windowpane, winter flounder, witch flounder and yellowtail flounder.

Code
# Name the list so it carries through
names(vonbert_species) <- vonbert_species

##### Reorganize species data into the list, and make group id for later
vonbert_species_agedat <-  map_dfr(vonbert_species, function(species_name){
  # Filter to the vonbert species
  regime_dat %>%
    filter(comname == species_name) %>% 
    mutate(group_id = str_c(regime, comname, sep = "-")) 
}, .id = "common name") 

# Species we want panels for:
key_species <- sort(c(
  "american plaice",
  "atlantic cod",
  "atlantic herring",
  "butterfish",
  "haddock", 
  "pollock",
  "scup",
  "silver hake",
  "summer flounder",
  "winter flounder",
  "yellowtail flounder"
))

Growth by Temperature Regime

To assess the impact of an elevated temperature regime, the size-at-age data was grouped into decadal periods relating to the shift in surface temperature regimes (2000-2009 & 2010-2019). Growth in fishes is commonly modeled using the “Von-Bertalanffy” growth function (VBGF) to capture how a fishes size (length or weight) changes with age. For each species the growth was modeled by fitting the VGBF to the length distribution-at-age.

\[L_t~=~L_{\infty}(1-e^{-K(t-t_0)})\]

Where \(L_t\) is the length in cm of a species at age \(t\). \(L_\infty\) is the asymptotic maximum length, \(K\) is _____, and \(t_0\) is x intercept, or the hypothetical age at 0 cm length.

Code
# Set the function for nls() to solve
vb <- vbFuns(param="Typical") 
Code
# Function to Pull Vonbert Coefficients, 
# In development: boot strapped intervals
vonbert_coef <- function(length_age_dat, start_points, min_obs = 20){
  
  ####  Von Bert Fitting Using: {FSA}  ####
  # Don't run on fewer than a minimum number of observations
  num_aged <- nrow(length_age_dat)
  
  # If number aged less than the minimum observations:
  # fill output with NA's, but preserve structure to build table
  if(num_aged < min_obs){
    
    # dataframe structure for coefficients
    na_coef_df <- data.frame(
      "Linf" = NA,
      "K" = NA,
      "t0" = NA)
    
    # dataframe structure for parameters
    na_param_df <- data.frame(
      "parameter" = NA,
      "estimate" = NA,
      "std_error" = NA,
      "t_value" = NA,
      "pr_t" = NA,
      "converged" = NA)
    
    # What to return when it fails/aborts:
    return(list("model"  = NA, 
                "coef"   = na_df, 
                "params" = na_param_df))
    
  } # Close conditional for sample size
  
  
  ####  Fitting Model When Sample Size is Sufficient  ####
  
  # Use nls to estimate VBGF parameters using starting points
  vbert_fit <- nls(length_cm ~ vb(age, Linf, K, t0),
                   data = length_age_dat,
                   start = start_points)

  # Access parameters with coef()
  vbert_coef <- as.data.frame(t(coef(vbert_fit)))
  
  
  # Get Summary Info (For convergence and std. error)
  vbert_summ <- summary(vbert_fit)
  
  # Convergence Info
  vbert_conv <- vbert_summ$convInfo$isConv
  
  # Parameters w/ Standard Error
  vbert_params <- as.data.frame(vbert_summ$parameters) %>% 
    rownames_to_column(var = "parameter") %>% 
    janitor::clean_names()
  
  # add convergence information in
  vbert_params$converged <- vbert_conv
  

  ####  Work in Progress / Debugging  ####
  # This section fails, can't find namespace/object "length_age_dat"
  # # Use bootstrapping for confidence intervals:
  # vbert_boot <- Boot(vbert_fit)  # Be patient! Be aware of some non-convergence
  # vbert_conf <- confint(vbert_boot)
  
  # output list
  out_l <- list(
    "model" = vbert_fit,
    "coef"  = vbert_coef,
    "params" = vbert_params#,
    #"conf"  = vbert_conf
  )

  # Return the list
  return(out_l)
  
}
Code
# Second Function for running that for a single Species, split by a factor column

# Function does these things:
# 1. grab data for one species
# 2. Get starting points for VB using all data on that species
# 3. Fit Vonbert to  both Regimes
process_group_vbert <- function(spec_name, split_col, min_obs){
  
  # Get the data using the species
  spec_dat <- dplyr::filter(vonbert_species_agedat, comname == spec_name)
  
  # Get starting points from all data
  test_starts <-  vbStarts(length_cm ~ age, 
                           data = spec_dat)
  
  # Get coefficients for groups
  spec_dat %>% 
    split(.[split_col]) %>% 
    map(.f = vonbert_coef,
            start_points = test_starts, 
            min_obs = min_obs) %>% 
    # map_dfr(pluck, "coef", .id = split_col) %>%  # Pull  out just the coefficients
    map_dfr(pluck, "params", .id = split_col) %>%  # Pull out the parameters and stderr
    mutate(comname = spec_name)
}
Code
# Running Each Species * Regime
species_coef <- vonbert_species %>% 
  map_dfr(function(x){
    process_group_vbert(
      spec_name = x, 
      split_col =  "regime", 
      min_obs = 30)
    
  })

Growth Impacts on Productivity

Yield-per-recruit analysis similar to Baudron et al. 2014?

Results

Temperature Regime Shift

Code
# Pull the whole area's temps
survey_temps <- filter(temp_regimes, survey_area == "Full Region")

# Pull the regime averages out
regime_avgs <- survey_temps %>% 
  distinct(regime, avg_temp_c, avg_anom_c, temp_label) 

# Get the regime averages separately for in text values
regime_1 <- regime_avgs %>% filter(regime == "1982-2009 Average")
regime_2 <- regime_avgs %>% filter(regime == "2010-2021 Average")


# Functions developed in 10_change_points.R
library(bcp)

# run bcp and flag changepoints using threshold
flag_bcp <- function(response_var,  
                     predictor_var = NULL, 
                     id_var, 
                     burnin = 10000,
                     mcmc = 10000, 
                     return.mcmc = T, 
                     threshold_prob = 0.3){
 
  # Use bcp to find support for change points
  bcp_x <- bcp(y = response_var, 
               x = predictor_var, 
               burnin = burnin, 
               mcmc = mcmc, 
               return.mcmc = return.mcmc)
  
  sink("/dev/null")
  # Convert summary to a dataframe
  bcp_sum <- as.data.frame(summary(bcp_x))
  sink()
  
  # Label the id variable
  bcp_sum$id <- id_var
  
  #Only pull records with support above desired threshold
  flag_dat <- bcp_sum[which(bcp_x$posterior.prob > threshold_prob), ]
  flag_dat <- as.data.frame(flag_dat)
  return(flag_dat)
  
}


# Set seed
set.seed(3905)


# Perform BCP regions temp breakpoints
anom_bcp <- flag_bcp(response_var = survey_temps$sst_anom, 
                     id_var = survey_temps$yr, 
                     threshold_prob = 0.3)

# in-text values
prob_text <- paste(round(anom_bcp$Probability *100, 2), collapse = ", ")

SST anomaly time series of the region indicate a jump in annual temperature anomalies indicative of a regime shift centered around 2010. A bayesian change point analysis showed highest support for a regime shift among 3 the years: 2009-2011. Change point probabilities for these years were 33.66, 32.51, 32.49% respectively.

Code
# Function to plot bcp breaks, highlighting their support
plot_bcp_breaks <- function(bcp_data, 
                            bcp_breaks, 
                            region_label, 
                            y_label = "sizeSpectra Exponent (b)",
                            y_col = "b"){
  
  # Set label elevation
  label_y <- mean(as.data.frame(bcp_data)[, y_col], na.rm = T) - 1
  
  # plot with breakpoint
  ggplot(bcp_data, aes_string("yr", y_col)) + 
    geom_line() + 
    geom_point() + 
    geom_vline(data = bcp_breaks, aes(xintercept = id), linetype = 2, color = "gray50") +
    ggrepel::geom_label_repel(data = bcp_breaks, 
                              aes(x = id, y = label_y, label = str_c(round(Probability, 2), "%"))) +
    labs(x = "Year", 
         y = y_label, 
         title = str_c(region_label, " - Bayesian Change Points"))
}





# Plot the bcp breaks for the Sea surface temperature anomalies
plot_bcp_breaks(bcp_data = survey_temps, 
                bcp_breaks = anom_bcp, 
                y_label = "Temperature Anomalies", 
                region_label = "Northeastern US Continental Shelf", 
                y_col = "sst_anom")

Average temperatures following a break in 2010, during the period of 2010-2021, were on average 1.09°C warmer than the than the average over the previous 28 years from 1982-2009. Temperatures within each regime were relatively stable, indicative of a jump in mean temperatures rather than rate change between them.

Code
# Plot how each one exhibits a break in the temperature across the region
temp_regimes %>% 
  filter(survey_area == "Full Region") %>% 
  ggplot( aes(yr, sst_anom)) +
  geom_col(aes(fill = regime),  size = 0.75, alpha = 0.4) +
  scale_fill_gmri() +
  # Using geomtextpath for label
  geom_textpath(
    aes(x = yr, 
        y = avg_anom_c, 
        label = temp_label, 
        linecolor = regime,
        colour = regime), 
        size = 4, vjust = -1, fontface = "bold", show.legend = FALSE) +
  # Line segments for the regime averages:
  geom_segment(
    aes(x = x, 
        xend = xend,
        y = avg_anom_c, 
        yend = avg_anom_c,
        color = regime),
    #key_glyph = "rect",
    size = .8,
    linetype = 1) +
  scale_x_continuous(expand = expansion(add = c(0.25,0.25))) +
  scale_y_continuous(labels = number_format(suffix = " \u00b0C")) +
  scale_color_gmri() +
  guides(fill = "none") +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.background = element_rect(fill = "transparent"),
    legend.key = element_rect(fill = "transparent", color = "transparent")) + 
  labs(
    title = "Northwest-Atlantic Temperature Regime Shift",
    x = "", 
    y = "Surface Temperature Anomaly",
    caption = "Anomalies calculated using 1982-2011 reference period.") +
  guides(label = "none",
         color = guide_legend(keyheight = unit(0.5, "cm"), ))

Distribution Shifts

Code
# Get the values for each of them in a standardized format
# Standardize to global benchmarks

# Latitude
lat_yrly <- biomass_weighted_var(station_totals, decdeg_beglat, est_year) %>% 
  left_join(select(global_benchmarks, c("comname", "avg_lat", "lat_sd")), by = "comname") %>% 
  mutate(biom_z = (biom_mu - avg_lat) / lat_sd) 

# Longitude
lon_yrly <- biomass_weighted_var(station_totals, decdeg_beglon, est_year) %>% 
  left_join(select(global_benchmarks, c("comname", "avg_lon", "lon_sd")), by = "comname") %>% 
  mutate(biom_z = (biom_mu - avg_lon) / lon_sd) 

# Bottom Temperature
t_yrly <- biomass_weighted_var(station_totals, bottemp, est_year) %>% 
  left_join(select(global_benchmarks, c("comname", "avg_temp", "temp_sd")), by = "comname") %>% 
  mutate(biom_z = (biom_mu - avg_temp) / temp_sd) 

# Depth
d_yrly <- biomass_weighted_var(station_totals, avgdepth, est_year) %>% 
  left_join(select(global_benchmarks, c("comname", "avg_depth", "depth_sd")), by = "comname") %>% 
  mutate(biom_z = (biom_mu - avg_depth) / depth_sd) 



# Put them together in one table for heatmap
standardized_changes <- list(
  "Latitude" = lat_yrly,
  "Longitude" = lon_yrly,
  "Bottom Temperature" = t_yrly,
  "Depth" = d_yrly) %>% 
  map_dfr(function(x){select(x, comname, est_year, var, biom_z)}, .id = "var_neat")



# Organize them by their global average latitude
global_benchmarks <- global_benchmarks %>% mutate(comname = fct_reorder(.f = comname, .x = avg_lat, .fun = median))

# Relevel the standardized changes
standardized_changes <- standardized_changes %>% 
  mutate(comname = factor(comname, 
                          levels = levels(global_benchmarks$comname)))


# Heatmap
standardized_changes %>% 
  ggplot(aes(est_year, comname, fill = biom_z)) +
  geom_tile() +
  scale_x_continuous(expand = expansion(add = c(0, 0))) +
  facet_wrap(~var_neat, ncol = 2, scales = "fixed") +
  scale_fill_distiller(palette = "RdBu",
                       oob = scales::oob_squish,
                       limits = c(-1.5, 1.5), 
                       breaks = c(-1.5, -.75, 0, .75, 1.5),
                       labels =c("<-1.5", "-.75", "0", ".75", ">1.5")) +
  labs(x = "Year", y = "",
       title = "Center of Biomass Movements Through Time",
       fill = "Center of Biomass Metric Change (z)",
       caption = "Species ordered based on center of biomass latitude across all years.") +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 7)) +
  guides(fill = guide_colorbar(title.position = "top", title.hjust = 0.5, barwidth = unit(6, "cm")))

Latitude

Code
# Standardized Plots
lat_yrly %>% 
  ggplot(aes(est_year, biom_z, group = comname)) +
  # geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +
  geom_smooth(method = "lm", formula = y ~ x, alpha = 0.8, se = FALSE, color = "gray50") +
  geom_jitter(width = 0.1, size = 2, alpha = 0.25) +
  geom_hline(aes(yintercept = 0, color = "Overall Average"), size = 1, linetype = 1) +
  scale_color_gmri(reverse = T) +
  labs(
    title = "Shift in Center of Biomass Latitudes",
    subtitle = "Yearly averages weighted by biomass, scaled against all-year average.",
    y = "Center of Biomass Latitude (z)",  
    x = "", 
    color = "Time Period", 
    fill = "Time Period") +
  theme(legend.position = "bottom")

Longitude

Code
# Longitude
lon_yrly %>% 
  ggplot(aes(est_year, biom_z, group = comname)) +
  # geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +
  geom_smooth(method = "lm", formula = y ~ x, alpha = 0.8, se = FALSE, color = "gray50") +
  geom_jitter(width = 0.1, size = 2, alpha = 0.25) +
  geom_hline(aes(yintercept = 0, color = "Overall Average"), size = 1, linetype = 1) +
  scale_color_gmri(reverse = T) +
  #coord_flip() +
  labs(
    title = "Shift in Center of Biomass Longtitudes",
    subtitle = "Yearly averages weighted by biomass, scaled against all-year average.",
    y = "Center of Biomass Longitude (z)",  
    x = "", 
    color = "Time Period", 
    fill = "Time Period") +
  theme(strip.text.y = element_text(size = 7, angle = 0, hjust = 0),
        legend.position = "bottom")

Mean Temperature

Code
# Temperature
t_yrly %>% 
  ggplot(aes(est_year, biom_z, group = comname)) +
  # geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +
  geom_smooth(method = "lm", formula = y ~ x, alpha = 0.8, se = FALSE, color = "gray50") +
  geom_jitter(width = 0.1, size = 2, alpha = 0.25) +
  geom_hline(aes(yintercept = 0, color = "Overall Average"), size = 1, linetype = 1) +
  scale_color_gmri(reverse = T) +
  #coord_flip() +
  labs(
    title = "Shift in Center of Biomass Temperature",
    subtitle = "Yearly averages weighted by biomass, scaled against all-year average.",
    y = "Center of Biomass Bottom Temperature (z)",  
    x = "", 
    color = "Time Period", 
    fill = "Time Period") +
  theme(legend.position = "bottom")

Mean Depth

Code
# Depth
d_yrly %>% 
  ggplot(aes(est_year, biom_z, group = comname)) +
  # geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +
  geom_smooth(method = "lm", formula = y ~ x, alpha = 0.8, se = FALSE, color = "gray50") +
  geom_jitter(width = 0.1, size = 2, alpha = 0.25) +
  geom_hline(aes(yintercept = 0, color = "Overall Average"), size = 1, linetype = 1) +
  scale_color_gmri(reverse = T) +
  #coord_flip() +
  labs(
    title = "Shift in Center of Biomass Depth",
    subtitle = "Yearly averages weighted by biomass, scaled against all-year average.",
    y = "Center of Biomass Depth (z)",  
    x = "", 
    color = "Time Period", 
    fill = "Time Period") +
  theme(strip.text.y = element_text(size = 7, angle = 0, hjust = 0),
        legend.position = "bottom")

Changes to Growth

Code
# Take the first part from the previous plot:
# takes each species
# for age 0 - max_age, solve what the size would be based on vb fit
regime_curve_prep <- function(spec_name){
  
  # The raw size/age data:
  t_dat <- vonbert_species_agedat %>% 
    filter(comname == spec_name)  
  
  # The group coefficients, augmented with x values across a range
  age_max <- max(t_dat$age)
  x_range <- seq(from = 0, to = age_max, by = .1)
  t_coef <- species_coef %>% 
    distinct(comname, regime, parameter, estimate) %>% 
    pivot_wider(values_from = "estimate", names_from = "parameter") %>% 
    filter(comname == spec_name) %>% 
    right_join(
      # Make a dataframe that covers the size range for the species
      y = data.frame(
        regime = rep(unique(t_dat$regime), length(x_range)),
        age = rep(x_range, length(unique(t_dat$regime)))),
      by = "regime") %>% 
    mutate(
      # Calculate the length at age
      pred_length = Linf * (1  - exp(-1 * K * (age - t0)))
    )
  
  # Pivot wider?
  # Take the size at age from each regime and get the difference
  t_coef <- t_coef %>%
    select(-c(Linf, K, t0)) %>% 
    pivot_wider(names_from = "regime", values_from = "pred_length") %>%
    mutate(
      growth_shift = `Warm Regime` - `Early Regime`,
      growth_frac = (growth_shift / `Early Regime`),
      size_change  = ifelse(growth_shift > 0, "Increase", "Decrease")
    )
  
  return(t_coef)
  
}


# Do it for all the species
regime_comparison <- vonbert_species %>% 
  map_dfr(.f = regime_curve_prep, .id = "comname")

Changes in the size-at-age growth relationship varied across species, with the majority (14/18) of the species exhibiting smaller L-infinity during the warmer temperature regime.

Code
# Try geom_waffle to just count each species as an individual and show that they are growing to smaller sizes in new regime:
library(ggwaffle)
age_set <- 4
ggplot(data = waffle_iron(
  data = filter(regime_comparison, age == age_set), 
  aes_d(group = size_change), rows = 2), 
         aes(x, y, fill = fct_rev(group))) +
  geom_waffle() +
  theme_waffle() +
  scale_fill_gmri() +
  theme(axis.text = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.title = element_text(face = "bold", size = 12),
        legend.position = "bottom",
        panel.grid.major.y = element_blank()) +
  labs(x = "", y = "", fill = str_c("Warm Regime Growth Impacts on ", age_set," Body Length:"))

The species that showed larger asymptotic lengths include the three hake species: red hake, silver hake, & white hake as well as Atlantic herring. These species are all among the majority group of species that is found primarily in the Northern part of the survey area, in and around the Gulf of Maine.

Code
# Reshape things so you can show confidence intervals here
plot_dat <- species_coef %>% 
  filter(parameter != "t0") %>% 
  mutate(llim = estimate - 2 * std_error,
         rlim = estimate + 2 * std_error,
         residence = ifelse(
           comname %in% c("atlantic mackerel", "scup", "black sea bass", "butterfish"),
           "Southern\nResident", 
           "Gulf of Maine\nResident"),
         regime = fct_rev(regime),
         comname = fct_rev(comname))


# Plot it
ggplot(plot_dat, aes(x = estimate, y = comname, group = regime, color = regime)) + 
    # need to nudge groups separately
  # Early Regime
    geom_segment(data = filter(plot_dat, regime == "Early Regime"), 
                 aes(x = llim, xend = rlim, yend = comname), 
                 position = position_nudge(x = 0, y = 0.25)) +
  geom_segment(data= filter(plot_dat, regime == "Warm Regime"), 
               aes(x = llim, xend = rlim, yend = comname), 
               position = position_nudge(x = 0, y = -0.25)) +
    geom_point(aes(color = regime), position = position_dodge(width = 1), size = 2) +
    #facet_grid(. ~ parameter, scales = "free", space = "free_y") +
    facet_grid(residence ~ parameter, scales = "free", space = "free_y") +
    scale_color_gmri(reverse = F) +
    labs(x = "Coefficient Estimate", y = "", fill = "Temperature Regime:") +
    theme(legend.position = "bottom", 
          panel.background = element_rect(color = "gray80", fill = "transparent", size = 1),
          panel.spacing = unit(2, "lines"),
          panel.border = element_rect(color = "gray80", fill = "transparent", size = 1))

Code
# Plotting Everything
species_coef %>%
  distinct(comname, regime, parameter, estimate) %>% 
  pivot_wider(values_from = c("estimate"), names_from = "parameter") %>% 
  pivot_longer(names_to = "vbcoef", 
               values_to = "val", 
               cols = c("K", "Linf", "t0")) %>% 
  filter(vbcoef != "t0") %>% 
  mutate(residence = ifelse(comname %in% c("atlantic mackerel", "scup", "black sea bass", "butterfish"), "Southern\nResident", "Gulf of Maine\nResident")) %>% 
  ggplot(aes(x = val, y = fct_rev(comname), fill = fct_rev(regime))) +
    geom_col(position = "dodge") +
    # facet_grid(. ~ vbcoef, scales = "free", space = "free_y") +
    facet_grid(residence ~ vbcoef, scales = "free", space = "free_y") +
    scale_fill_gmri(reverse = T) +
    labs(x = "Coefficient Estimate", y = "", fill = "Temperature Regime:") +
  theme(legend.position = "bottom", 
        panel.background = element_rect(color = "gray80", fill = "transparent", size = 1),
        panel.spacing = unit(2, "lines"),
        panel.border = element_rect(color = "gray80", fill = "transparent", size = 1))

Changes in \(L_\infty\)

Code
species_coef %>% 
  distinct(comname, regime, parameter, estimate) %>% 
  pivot_wider(values_from = "estimate", names_from = "parameter") %>% 
  mutate(regime_2 = ifelse(regime == "Early Regime", "early", "warm")) %>% 
  select(comname, regime_2, Linf) %>% 
  pivot_wider(values_from = Linf, names_from = regime_2) %>% 
  mutate(linf_diff = warm - early,
         diff_perc = (linf_diff / early) * 100) %>% 
  mutate(across(where(is.numeric), round, 2)) %>% 
  arrange(desc(diff_perc)) %>% 
  setNames(c("Common Name", "Early Regime Linf", "Warm Regime Linf", "Change in Linf", "Percent Change")) %>% 
  gt()
Common Name Early Regime Linf Warm Regime Linf Change in Linf Percent Change
silver hake 49.85 59.58 9.73 19.51
white hake 164.88 187.04 22.16 13.44
red hake 43.67 46.28 2.61 5.97
windowpane 30.05 31.36 1.31 4.36
atlantic herring 29.13 29.56 0.44 1.50
atlantic mackerel 45.60 45.77 0.17 0.37
acadian redfish 36.94 36.53 -0.41 -1.12
black sea bass 78.72 77.42 -1.30 -1.65
pollock 128.71 123.33 -5.38 -4.18
butterfish 26.30 24.80 -1.50 -5.71
haddock 64.14 59.66 -4.48 -6.98
summer flounder 84.12 77.44 -6.68 -7.95
yellowtail flounder 47.77 43.55 -4.22 -8.84
winter flounder 52.17 47.18 -4.99 -9.57
atlantic cod 122.77 110.31 -12.46 -10.15
witch flounder 51.76 45.96 -5.81 -11.22
american plaice 57.39 44.40 -12.98 -22.63
scup 68.32 42.36 -25.96 -38.00

VB Fits to Data:

Code
# Make the plot for each species
vb_fits_plot = function(spec_name){
    
  # The raw size/age data:
  t_dat <- vonbert_species_agedat %>% filter(comname == spec_name)  
  
  # The group coefficients, augmented with x values across a range
  age_max <- max(t_dat$age)
  x_range <- seq(from = 0, to = age_max, by = .1)
  t_coef <- species_coef %>% 
    distinct(comname, regime, parameter, estimate) %>% 
    pivot_wider(values_from = "estimate", names_from = "parameter") %>% 
    filter(comname == spec_name) %>% 
    right_join(
      y = data.frame(
        regime = rep(unique(t_dat$regime), length(x_range)),
        age = rep(x_range, length(unique(t_dat$regime)))),
      by = "regime") %>% 
    mutate(
      # Calculate the length at age
      pred_length = Linf * (1  - exp(-1 * K * (age - t0)))
    )
  
  
  # Make plot?
  ggplot(t_dat, aes(x = age, y = length_cm)) +
    geom_point(alpha = 0.25) +
    #facet_wrap(~regime, ncol = 1) +
    geom_line(data = t_coef,
              aes(x = age, y = pred_length, group = regime, color = regime), 
              size = 1) +
    scale_x_continuous(breaks = seq(0, age_max, by = 2),
                       limits = c(0, age_max),
                       expand = expansion(mult = c(0.1, .2))) +
    scale_color_gmri() +
    geom_dl(data = t_coef, 
            aes(x = age, y = pred_length, label = regime, color = regime), 
            method= list("last.qp", cex = 0.65, hjust = -.15),
            size = 12) +
    labs(x = "Age",
         y = "Length (cm)",
         title = spec_name,
         color = "Temperature Regime:",
         subtitle = "Length at Age Regimes")

    
  }
Code
# Loop through the key species:
walk(key_species, 
     plot_panelset, 
     plot_fun = vb_fits_plot)

american plaice

atlantic cod

atlantic herring

butterfish

haddock

pollock

scup

silver hake

summer flounder

winter flounder

yellowtail flounder

Discussion


Bonus: Data Explorations

Age Data Availability

The availability of age data varied by species and through time, with older ages becoming less frequently sampled in recent years across many species.

Code
# Plotting age distribution as ggridges
plot_age_ridgeplot <- function(species){
  vonbert_species_agedat %>% 
    filter(comname == species,
           age <= 15) %>% 
    mutate(yr = factor(est_year),
           yr = fct_rev(yr)) %>% 
    ggplot(aes(x = age, y = yr)) +
      geom_density_ridges(aes(fill = regime), color = "white", show.legend = FALSE) +
      scale_fill_gmri() +
      facet_wrap(~comname) +
      scale_x_continuous(limits = c(0,NA), 
                         expand = expansion(add = c(0, 0))) +
      labs(x = "Age", y = "Year") +
      guides(size = guide_legend(title.position = "top", title.hjust = 0.5))
  
}
Code
# Loop through the key species:
walk(key_species, plot_panelset, plot_fun = plot_age_ridgeplot)

american plaice

atlantic cod

atlantic herring

butterfish

haddock

pollock

scup

silver hake

summer flounder

winter flounder

yellowtail flounder

Strong age classes were particularly prominent in the following species, showing lasting numbers in subsequent years.

Species Strong Year Classes
Atlantic Cod 2003, 2013
Haddock 2003, 2010, 2013
Atlantic Herring 2008, 2011
Pollock 2000, 2012

Growth Increment Changes

Annual Growth increments for each species were computed to track the change in growth of specific age-groups through time.

Code
# Age-Size Increments Differences
vb_incs <- vonbert_species_agedat %>% 
  group_by(`common name`, est_year,  age, regime) %>% 
  summarise(avg_len = mean(length_cm, na.rm = T),
            avg_wt  = mean(indwt, na.rm = T),
            .groups = "drop") %>% 
  arrange(`common name`, age) %>% 
  group_split(`common name`, est_year) %>% 
  map_dfr(function(x){
    x %>% 
      mutate(
        # What age is being subtracted from
        lead_age   = lead(age, n = 1),
        # String built from the value to confirm
        diff_name  = str_c(age, " - ", lead_age),
        # Length/Weight at age + 1
        len_age_p1 = lead(avg_len, n = 1),
        wt_age_p1  = lead(avg_wt, n = 1),
        # Length/Weight Growth Increments
        len_inc    = len_age_p1 - avg_len,
        wt_inc     = wt_age_p1 - avg_wt,
        # Percent change in bodymass for going from age to age+1
        len_perc_change = (len_inc / avg_len) * 100,
        wt_perc_change = (wt_inc / avg_wt) * 100) %>% 
      filter((lead_age - age == 1))
  })

Age specific growth (annual growth increments) help detail where in a species life cycle it may be experiencing favorable or less-favorable conditions for growth. These can be a result of differences in habitat use or prey resources at different ages, and/or a change in exposure to temperature or other stressors that inflict a metabolic cost.

Code
# Put together a timeseries of:
# how big were the age 0's,
# how big were the age 1's

# Plotting Age Increments
plot_age_increments <- function(
    spec_choice, 
    var_choice = "length", 
    gam_knots = 15, 
    facet_prefix = "Age Increment: ",
    facet_suffix = ""){
  
# Code to switch variables and their labels:
  var_sym <- switch(var_choice,
    "length" = sym("len_inc"),
    "weight" = sym("wt_inc"))
  var_label <- switch(var_choice,
                      "length" = "Annual Growth  (cm)",
                      "weight" = "Annual Growth  (kg)")
  
  
  # Add pre/suffix to facet labels
  appender <- function(string, prefix = facet_prefix, suffix = facet_suffix) {
    paste0(prefix, string)}

  # Make the plot
  vb_incs %>% 
    filter(`common name` == spec_choice,
           lead_age <= 5) %>% 
    drop_na({{var_sym}}) %>% 
    ggplot(aes(est_year, {{var_sym}})) +
      # geom_col(fill = gmri_cols("gmri blue")) +
      geom_point(color = gmri_cols("gmri blue")) +
      geom_line(color = gmri_cols("gmri blue"), group = 1, alpha = 0.3, size = 0.8) +
      geom_smooth(method = "gam",
                  formula = y ~ s(x, bs = "cs", k = gam_knots),
                  se = FALSE,
                  size = 1,
                  color = gmri_cols("orange")) +
      facet_wrap(~diff_name, ncol = 1, scales = "free",
                 labeller = as_labeller(appender))  +
      theme(legend.position = "bottom") +
      labs(x = "",
           y = var_label,
           color = "Age Increment",
           title = str_to_title(spec_choice),
           subtitle = str_c("Changes in annual ", var_choice, " growth increments."))

}
Code
# Loop through the key species:
walk(key_species, plot_panelset, plot_fun = plot_age_increments, var_choice = "weight", gam_knots = 5)

american plaice

atlantic cod

atlantic herring

butterfish

haddock

pollock

scup

silver hake

summer flounder

winter flounder

yellowtail flounder

Length at Age Changes

Code
# Distribution of Length/Weight Across Regimes 
plot_regime_sizes <- function(spec, 
                              var = "length", 
                              age_lim = 5, 
                              facet_prefix = "Age - ", 
                              facet_suffix = ""){
  
  # switches for length or width
  var_col <- switch (var,
    "length" = sym("length_cm"),
    "weight" = sym("indwt"))
  var_lab <- switch (var,
    "length" = "Individual Length (cm)",
    "weight" = "Individual Weight (kg)")
  
  # Add pre/suffix to facet labels
  appender <- function(string, prefix = facet_prefix, suffix = facet_suffix) {
    paste0(prefix, string)}
  
  # Select Ages, Prep data
  size_changes <- regime_dat %>% 
  filter(
    comname == spec,
    age <= 6) %>% 
  ## Calculate group averages
  group_by(age, regime) %>% 
  mutate(
    median_size = median(!!{{var_col}})) %>% 
    ungroup() 
  
  # Make median dataframe
  median_df <- size_changes %>%
    distinct(age, regime, median_size) %>%
    pivot_wider(names_from = "regime", values_from = "median_size")
  
  # Make Plot
  size_changes %>% 
    ggplot(aes(x = !!{{var_col}})) +
    
    ## distribution
    stat_halfeye(aes(color = regime, 
                     fill = after_scale(color)),
                 slab_alpha = .5, 
                 .width = 0, 
                 trim = T, 
                 shape = 21, 
                 point_colour = "gray25",
                 stroke = 1.6,
                 scale = .86) +
    
    # median line
    geom_segment(
      data = median_df,
      aes(x = `Early Regime`,
          xend = `Warm Regime`,
          y = 0,
          yend = 0),
      size = .8,
      color = "grey25", # , position = position_nudge(y = -.15)
    ) +
    
    ## median points - not working
    stat_halfeye(aes(color = regime), .width = c(0), slab_fill = NA) +
    
   
    facet_wrap(~age, ncol = 1, 
               labeller = as_labeller(appender), scales = "free") +
    scale_color_gmri() + 
    scale_fill_gmri() +
    scale_y_continuous(labels = percent_format()) +
    theme(legend.position = "bottom") +
    labs(x = var_lab,
         y = "Proportion",
         title = spec,
         fill = "Temperature Regime",
         color = "Temperature Regime",
         shape = "Median"
         )
  
  
}


# Tester: 
# plot_regime_sizes(key_species[[1]], var = "length")
Code
# Loop through the key species:
walk(key_species, plot_panelset, plot_fun = plot_regime_sizes, var = "length")

american plaice

atlantic cod

atlantic herring

butterfish

haddock

pollock

scup

silver hake

summer flounder

winter flounder

yellowtail flounder

Weight at Age Changes

Code
# Loop through the key species:
walk(key_species, plot_panelset, plot_fun = plot_regime_sizes, var = "weight")

american plaice

atlantic cod

atlantic herring

butterfish

haddock

pollock

scup

silver hake

summer flounder

winter flounder

yellowtail flounder

Shifts in Growth

Code
# Plot their growth curves
regime_comparison %>% 
  #filter(comname %in% key_species) %>% 
  ggplot(aes(x = age, y = growth_frac)) +
    geom_hline(yintercept = 0, size = 0.5, lty = 2, color = "black") +
    geom_line(aes(group = comname), show.legend = FALSE, size = 1, color = "orange")+
    scale_y_continuous(labels = percent_format()) +
    facet_wrap(~comname, scales = "free_x", ncol = 4) + 
    labs(x = "Age", y = "Shift in Length Moving to Warm Regime", 
         title = "Thermal Regime Impacts on Growth",
         subtitle = "Changes in Size at Age (Warm Regime - Historic Regime)") 

Code
# Kathy requested the curve data.
# Going to save both the coefficients and the fitted data for all species
library(here)
write_csv(regime_comparison, here("data/growth_results/growth_regime_fits.csv"))
write_csv(species_coef, here("data/growth_results/growth_regime_coefs.csv"))
Code
# Try horizons
library(ggHoriPlot)


# Filter outliers beyond 5th and 95th quantile
cutpoints <- regime_comparison %>% 
  mutate(
    outlier = between(
      growth_frac,
      quantile(growth_frac, 0.05, na.rm = T)-
        1.5 * IQR(growth_frac, na.rm = T),
      quantile(growth_frac, 0.95, na.rm = T)+
        1.5 * IQR(growth_frac, na.rm = T))) %>% 
  filter(outlier)


# Set origin
ori <- 0


# Manually pick cut points
sca <- seq(-1, 1, length.out = 9)
sca_bins <- str_c(sca[1:(length(sca)-1)], " to ", sca[2:length(sca)], "%")
sca_labels <- rev(sca_bins)


theme_set(theme_bw() + 
            theme(
              # Titles
              plot.title = element_text(hjust = 0, face = "bold", size = 14),
              plot.subtitle = element_text(size = 9),
              plot.caption = element_text(size = 8, margin = margin(t = 20), color = "gray40"),
              plot.margin = unit(c(1, 1, 2, 1), "lines"),
              legend.title = element_text(size = 9),
              legend.text = element_text(size = 9),
              # Axes
              axis.line.y = element_line(color = "black"),
              axis.ticks.y = element_line(), 
              axis.line.x = element_line(color = "black"),
              axis.ticks.x = element_line(), 
              axis.text = element_text(size = 11),
              axis.title = element_text(size = 12),
              rect = element_rect(fill = "transparent", color = "black"),
              # Facets
              strip.text = element_text(color = "white", 
                                        face = "bold",
                                        size = 11),
              strip.background = element_rect(
                color = "#00736D", 
                fill = "#00736D", 
                size = 1, 
                linetype="solid")))


# Plot it all
ggplot(regime_comparison, aes(x = age, y = growth_frac)) +
  geom_horizon(
    aes(age, 
        y = growth_frac,
        fill = ..Cutpoints..), 
    origin = ori, 
    horizonscale = sca) +
  scale_fill_hcl(palette = 'RdBu', reverse = T, labels = sca_labels) +
  facet_grid(comname~.) +
  theme_bw() + 
  scale_x_continuous(limits = c(0,10), expand = expansion(add = c(0, 0))) +
  theme(
    # Titles
    plot.title = element_text(hjust = 0, face = "bold", size = 14),
    plot.subtitle = element_text(size = 9),
    plot.caption = element_text(
      size = 8, 
      margin = margin(t = 20), 
      color = "gray40"),
    plot.margin = unit(c(1, 1, 2, 1), "lines"),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 9),
    # Axes
    axis.line.y = element_line(color = "black"),
    axis.line.x = element_line(color = "black"),
    axis.ticks.x = element_line(), 
    axis.text = element_text(size = 11),
    axis.title = element_text(size = 12),
    rect = element_rect(fill = "transparent", color = "black"),
    # Facets
    strip.text = element_text(color = "white", 
                              face = "bold",
                              size = 11),
    strip.background = element_rect(
      color = "#00736D", 
      fill = "#00736D", 
      size = 1, 
      linetype="solid"),
    panel.grid = element_blank(),
    panel.spacing.y=unit(0.05, "lines"),
    strip.text.y = element_text(size = 7, angle = 0, hjust = 0, face = "bold"),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.border = element_blank(),
    legend.position = "left")  +
  labs(x = 'Age', 
       fill = "Change in Body Size at Age") +
  guides(fill = guide_legend(ncol = 1))

 

A work by Adam A. Kemberling

Akemberling@gmri.org